###################################################################
### Political Institutions, Energy Transitions, and Air Quality ### 
### Evidence from Global Urban Areas                            ###
### Main empirical findings / difference in differences         ###                  
### This version: November 16, 2024                             ###                   
###################################################################

#The code below reproduces the main empirical findings in the paper and SI related to the differences in differences, specifically

####################
###### Figures #####
####################

#Main Text

#Figure 13: Event study analysis for CO2 after the retirement of a coal power unit.
#Figure 14: Event Study Analysis for PM2.5 and NO2 after the Retirement of a Coal Power Plant Unit.
#Figure 15: Event study analysis for mortality rates associated with PM2.5.
#Figure 16: Event study analysis for mortality rates associated with NO2.
#Figure 17: Heterogeneous effects of phasing out coal power plants on public health outcomes (Mortality rates-particulate matter).
#Figure 18: Heterogeneous effects of phasing out coal power plants on public health outcomes (Mortality rates-nitrogen dioxide).

#Supplemental Information

#Figure S1: Staggered difference-in-differences for PM2.5 Rates
#Figure S2: Staggered difference-in-differences for NO2 Rates.
#Figure S3: Staggered differences-in-differences, cohort aggregation (PM2.5 Rates).
#Figure S4: Staggered differences-in-differences, cohort aggregation (PM2.5 Rates).
#Figure S5: Staggered differences-in-differences, cohort aggregation (PM2.5 Rates).
#Figure S6: Staggered differences-in-differences, cohort aggregation (NO2 Rates).
#Figure S7: Staggered differences-in-differences, cohort aggregation (NO2 Rates).
#Figure S8: Staggered differences-in-differences, cohort aggregation (NO2 Rates).
#Figure S9: Staggered differences-in-differences, event study (PM2.5 Rates).
#Figure S10: Staggered differences-in-differences, event study (PM2.5 Rates) by capacity.
#Figure S11: Staggered differences-in-differences, event study (NO2 Rates).
#Figure S12: Staggered differences-in-differences, event study (NO2 Rates) by capacity.
#Figure S13: FEct for Mortality Rates from PM2.5.
#Figure S14: FEct for Emissions of PM2.5.
#Figure S15: FEct for Mortality Rates from NO2.
#Figure S16: FEct for Emissions of NO2.
#Figure S17: FEct for Mortality Rates from PM2.5.
#Figure S18: FEct for Emissions of PM2.5.
#Figure S19: FEct for Mortality Rates from NO2.
#Figure S20: FEct for Emissions of NO2.
#Figure S21: Sensitivity Analysis PM2.5 Mortality
#Figure S22: Sensitivity Analysis PM2.5 Emissions
#Figure S23: Sensitivity Analysis NO2 Mortality
#Figure S24: Sensitivity Analysis NO2 Emissions

####################
###### Tables  #####
####################

#Main Text

#Table 1: Democracy and Air Pollution (Particulate Matter)
#Table 2: Democracy and Air Pollution (NO2)

#Supplemental Information

#Table S2: Standard TWFE: Mortality Rates.
#Table S3: Standard TWFE: Total Emissions.
#Table S4: Standard TWFE: Mortality Rates.
#Table S5: Standard TWFE: Total Emissions
#Table S6: Standard Two-Way Fixed Effects: Mortality Rates by Distance to the Closest Retired Unit
#Table S7: Standard Two-Way Fixed Effects: Mortality Rates by Distance to the Closest Retired Unit.
#Table S8: Standard Two-Way Fixed Effects: Emissions by Distance to the Closest Retired Unit.
#Table S9: Standard Two-Way Fixed Effects: Emissions by Distance to the Closest Retired Unit.
#Table S10: Standard Two-Way Fixed Effects: Emissions by Distance to the Closest Retired Unit.
#Table S11: Standard Two-Way Fixed Effects: Mortality Rates by Capacity of the Closest Retired Coal Unit.
#Table S12: Standard Two-Way Fixed Effects: Mortality Rates by Capacity of the Closest Retired Coal Unit. 
#Table S13: Standard Two-Way Fixed Effects: Emissions by Capacity of the Closest Retired Coal Unit.
#Table S14: Standard Two-Way Fixed Effects: Emissions by Capacity of the Closest Retired Coal Unit.
#Table S15: Standard Two-Way Fixed Effects: Emissions by Capacity of the Closest Retired Coal Unit
#Table S16: Standard Two-Way Fixed Effects: Full Sample of All Cities.
#Table S17: Standard Two-Way Fixed Effects: Full Sample of All Cities.

#---------------------------------------------------#
# Part 1: Preliminaries                             #
#---------------------------------------------------#

rm(list = ls())

#List of packages for session
.packages = c("tidyverse", "did", "texreg", "countrycode",
              "fixest", "WDI", "gridExtra", "sf", "nngeo", "fect")

#Install CRAN packages (if not already installed)
.inst = .packages %in% installed.packages()
if(length(.packages[!.inst]) > 0) install.packages(.packages[!.inst])

#Loading packages into session 
lapply(.packages, require, character.only=TRUE)

# Set Working Directory to wherever files are downloaded
#setwd()

#---------------------------------------------------#
# Part 2: Loading the data                          #
#---------------------------------------------------#

master_cities_panel <- read.csv("master_cities_panel.csv")
master_pollution_dem <- read.csv("master_pollution_dem.csv")
master_cities_panel_robustness <- read.csv("master_cities_panel_robustness.csv")
master_plants_panel <- read.csv("master_plants_panel.csv")

#---------------------------------------------------#
# Part 3: Preparing the data for the staggered did  #
#---------------------------------------------------#

#Creating data for PM2.5 Rates

did_master_cities_panel_pm <- master_cities_panel %>%
  mutate(Retired.year50k = ifelse(is.na(Retired.year50k), 0, Retired.year50k),
         Retired.year100k = ifelse(is.na(Retired.year100k), 0, Retired.year100k),
         Retired.year150k = ifelse(is.na(Retired.year150k), 0, Retired.year150k),
         Retired.year450k = ifelse(is.na(Retired.year450k), 0, Retired.year450k),
         Retired.year50k_high = ifelse(is.na(Retired.year50k_high), 0, Retired.year50k_high),
         Retired.year50k_low = ifelse(is.na(Retired.year50k_low), 0, Retired.year50k_low)) %>%
  dplyr::select(ID, year, Rates_PM, Retired.year50k, Retired.year100k, 
                Retired.year150k, Retired.year450k, Population, CO2, political_regime,
                Retired.year50k_high, Retired.year50k_low) %>% na.omit()

#Creating data for NO2 Rates

did_master_cities_panel_no2 <- master_cities_panel %>%
  mutate(Retired.year50k = ifelse(is.na(Retired.year50k), 0, Retired.year50k),
         Retired.year100k = ifelse(is.na(Retired.year100k), 0, Retired.year100k),
         Retired.year150k = ifelse(is.na(Retired.year150k), 0, Retired.year150k),
         Retired.year450k = ifelse(is.na(Retired.year450k), 0, Retired.year450k),
         Retired.year50k_high = ifelse(is.na(Retired.year50k_high), 0, Retired.year50k_high),
         Retired.year50k_low = ifelse(is.na(Retired.year50k_low), 0, Retired.year50k_low)) %>%
  dplyr::select(ID, year, Rates_NO2, Retired.year50k, Retired.year100k, 
                Retired.year150k, Retired.year450k, Population, CO2, political_regime,
                Retired.year50k_high, Retired.year50k_low) %>% na.omit()

#Creating data for PM Emissions

did_master_cities_panel_pm_emissions <- master_cities_panel %>%
  mutate(Retired.year50k = ifelse(is.na(Retired.year50k), 0, Retired.year50k),
         Retired.year100k = ifelse(is.na(Retired.year100k), 0, Retired.year100k),
         Retired.year150k = ifelse(is.na(Retired.year150k), 0, Retired.year150k),
         Retired.year450k = ifelse(is.na(Retired.year450k), 0, Retired.year450k),
         Retired.year50k_high = ifelse(is.na(Retired.year50k_high), 0, Retired.year50k_high),
         Retired.year50k_low = ifelse(is.na(Retired.year50k_low), 0, Retired.year50k_low)) %>%
  dplyr::select(ID, year, PM, Retired.year50k, Retired.year100k, 
                Retired.year150k, Retired.year450k, Population, CO2, political_regime,
                Retired.year50k_high, Retired.year50k_low) %>% na.omit()

#Creating data for NO2 Emissions

did_master_cities_panel_no2_emissions <- master_cities_panel %>%
  mutate(Retired.year50k = ifelse(is.na(Retired.year50k), 0, Retired.year50k),
         Retired.year100k = ifelse(is.na(Retired.year100k), 0, Retired.year100k),
         Retired.year150k = ifelse(is.na(Retired.year150k), 0, Retired.year150k),
         Retired.year450k = ifelse(is.na(Retired.year450k), 0, Retired.year450k),
         Retired.year50k_high = ifelse(is.na(Retired.year50k_high), 0, Retired.year50k_high),
         Retired.year50k_low = ifelse(is.na(Retired.year50k_low), 0, Retired.year50k_low)) %>%
  dplyr::select(ID, year, NO2, Retired.year50k, Retired.year100k, 
                Retired.year150k, Retired.year450k, Population, CO2, political_regime,
                Retired.year50k_high, Retired.year50k_low) %>% na.omit()

#Creating data for CO2 Emissions

did_master_cities_panel_co2_emissions <- master_cities_panel %>%
  mutate(Retired.year50k = ifelse(is.na(Retired.year50k), 0, Retired.year50k),
         Retired.year100k = ifelse(is.na(Retired.year100k), 0, Retired.year100k),
         Retired.year150k = ifelse(is.na(Retired.year150k), 0, Retired.year150k),
         Retired.year450k = ifelse(is.na(Retired.year450k), 0, Retired.year450k),
         Retired.year50k_high = ifelse(is.na(Retired.year50k_high), 0, Retired.year50k_high),
         Retired.year50k_low = ifelse(is.na(Retired.year50k_low), 0, Retired.year50k_low)) %>%
  dplyr::select(ID, year, NO2, Retired.year50k, Retired.year100k, 
                Retired.year150k, Retired.year450k, Population, CO2, political_regime,
                Retired.year50k_high, Retired.year50k_low) %>% na.omit() %>% filter(CO2 > 0) %>%
  mutate(logco2 = log(CO2))

###############################################
######### FIGURES MAIN TEXT ANALYSIS ##########
###############################################

#---------------------------------------------------#
# Figure 13: Event study analysis for CO2           #
# after the retirement of a coal power unit.        #
#---------------------------------------------------#

attgt_retirement50k_high_var <- att_gt(yname = "logco2",
                                       tname = "year",
                                       idname = "ID",
                                       gname = "Retired.year50k_high",
                                       base_period = "varying",
                                       bstrap = T,
                                       clustervars = "ID",
                                       control_group = "nevertreated",
                                       #xformla = ~Population+CO2,
                                       data = did_master_cities_panel_co2_emissions
)
simple.attgt_retirement50k_high_var <- aggte(attgt_retirement50k_high_var, type = "simple")
dynamic.attgt_retirement50k_high_var <- aggte(attgt_retirement50k_high_var, type = "dynamic")

attgt_retirement50k_low_var <- att_gt(yname = "logco2",
                                      tname = "year",
                                      idname = "ID",
                                      gname = "Retired.year50k_low",
                                      base_period = "varying",
                                      bstrap = T,
                                      clustervars = "ID",
                                      control_group = "nevertreated",
                                      #xformla = ~Population+CO2,
                                      data = did_master_cities_panel_co2_emissions
)
simple.attgt_retirement50k_low_var <- aggte(attgt_retirement50k_low_var, type = "simple")
dynamic.attgt_retirement50k_low_var <- aggte(attgt_retirement50k_low_var, type = "dynamic")

mod1_pm_retirement50k_high_var <- tidy(dynamic.attgt_retirement50k_high_var) %>% mutate(model = "50k from Retired Coal Unit (High Capacity MW)")
mod1_pm_retirement50k_low_var <- tidy(dynamic.attgt_retirement50k_low_var) %>% mutate(model = "50k from Retired Coal Unit (Low Capacity MW)")

mod1_pm_retirement_robustness2 <- bind_rows(mod1_pm_retirement50k_high_var, mod1_pm_retirement50k_low_var)

pdf(file = "results/fig14c_co2emissions_capacity.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot(mod1_pm_retirement_robustness2, aes(x = as.factor(event.time), y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange(aes(color = model, group = model), position = position_dodge(width=0.3)) +
  xlab("Time to Coal Power Unit Retirement") +
  ylab("Coefficient") +
  ggtitle("Average Effect by Length of Exposure") +
  theme_classic() + scale_color_brewer(palette = "Set1")  + 
  geom_vline(xintercept = "0", linetype = "twodash", color = "grey50", size=0.5) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", size=0.3)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure 14: Event Study Analysis for PM2.5 and NO2 after the Retirement of a Coal Power Plant Unit                  #
#--------------------------------------------------------------------------------------------------------------------#

#Left-Panel

attgt_retirement50k_high_var_pm <- att_gt(yname = "PM",
                                       tname = "year",
                                       idname = "ID",
                                       gname = "Retired.year50k_high",
                                       base_period = "varying",
                                       bstrap = T,
                                       clustervars = "ID",
                                       control_group = "nevertreated",
                                       #xformla = ~Population+CO2,
                                       data = did_master_cities_panel_pm_emissions
)
simple.attgt_retirement50k_high_var_pm <- aggte(attgt_retirement50k_high_var_pm, type = "simple")
dynamic.attgt_retirement50k_high_var_pm <- aggte(attgt_retirement50k_high_var_pm, type = "dynamic")

attgt_retirement50k_low_var_pm <- att_gt(yname = "PM",
                                      tname = "year",
                                      idname = "ID",
                                      gname = "Retired.year50k_low",
                                      base_period = "varying",
                                      bstrap = T,
                                      clustervars = "ID",
                                      control_group = "nevertreated",
                                      #xformla = ~Population+CO2,
                                      data = did_master_cities_panel_pm_emissions
)
simple.attgt_retirement50k_low_var_pm <- aggte(attgt_retirement50k_low_var_pm, type = "simple")
dynamic.attgt_retirement50k_low_var_pm <- aggte(attgt_retirement50k_low_var_pm, type = "dynamic")

mod1_pm_retirement50k_high_var <- tidy(dynamic.attgt_retirement50k_high_var_pm) %>% mutate(model = "50k from Retired Coal Unit (High Capacity MW)")
mod1_pm_retirement50k_low_var <- tidy(dynamic.attgt_retirement50k_low_var_pm) %>% mutate(model = "50k from Retired Coal Unit (Low Capacity MW)")

mod1_pm_retirement_robustness2 <- bind_rows(mod1_pm_retirement50k_high_var, mod1_pm_retirement50k_low_var)

pdf(file = "results/fig12c_pmemissions_capacity.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot(mod1_pm_retirement_robustness2, aes(x = as.factor(event.time), y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange(aes(color = model, group = model), position = position_dodge(width=0.3)) +
  xlab("Time to Coal Power Unit Retirement") +
  ylab("Coefficient") +
  ggtitle("Average Effect by Length of Exposure") +
  theme_classic() + scale_color_brewer(palette = "Set1")  + 
  geom_vline(xintercept = "0", linetype = "twodash", color = "grey50", size=0.5) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", size=0.3)
dev.off()

#Right-Panel

attgt_retirement50k_high_var_no2 <- att_gt(yname = "NO2",
                                       tname = "year",
                                       idname = "ID",
                                       gname = "Retired.year50k_high",
                                       base_period = "varying",
                                       bstrap = T,
                                       clustervars = "ID",
                                       control_group = "nevertreated",
                                       #xformla = ~Population+CO2,
                                       data = did_master_cities_panel_no2_emissions
)
simple.attgt_retirement50k_high_var_no2 <- aggte(attgt_retirement50k_high_var_no2, type = "simple")
dynamic.attgt_retirement50k_high_var_no2 <- aggte(attgt_retirement50k_high_var_no2, type = "dynamic")

attgt_retirement50k_low_var_no2 <- att_gt(yname = "NO2",
                                      tname = "year",
                                      idname = "ID",
                                      gname = "Retired.year50k_low",
                                      base_period = "varying",
                                      bstrap = T,
                                      clustervars = "ID",
                                      control_group = "nevertreated",
                                      #xformla = ~Population+CO2,
                                      data = did_master_cities_panel_no2_emissions
)
simple.attgt_retirement50k_low_var_no2 <- aggte(attgt_retirement50k_low_var_no2, type = "simple")
dynamic.attgt_retirement50k_low_var_no2 <- aggte(attgt_retirement50k_low_var_no2, type = "dynamic")

mod1_no2_retirement50k_high_var <- tidy(dynamic.attgt_retirement50k_high_var_no2) %>% mutate(model = "50k from Retired Coal Unit (High Capacity MW)")
mod1_no2_retirement50k_low_var <- tidy(dynamic.attgt_retirement50k_low_var_no2) %>% mutate(model = "50k from Retired Coal Unit (Low Capacity MW)")

mod1_no2_retirement_robustness2 <- bind_rows(mod1_no2_retirement50k_high_var, mod1_no2_retirement50k_low_var)

pdf(file = "results/fig13c_no2emissions_capacity.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot(mod1_no2_retirement_robustness2, aes(x = as.factor(event.time), y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange(aes(color = model, group = model), position = position_dodge(width=0.3)) +
  xlab("Time to Coal Power Unit Retirement") +
  ylab("Coefficient") +
  ggtitle("Average Effect by Length of Exposure") +
  theme_classic() + scale_color_brewer(palette = "Set1")  + 
  geom_vline(xintercept = "0", linetype = "twodash", color = "grey50", size=0.5) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", size=0.3)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure 15: Event study analysis for mortality rates associated with PM2.5                                          #
#--------------------------------------------------------------------------------------------------------------------#

attgt_retirement50k_var_pm <- att_gt(yname = "Rates_PM",
                                  tname = "year",
                                  idname = "ID",
                                  gname = "Retired.year50k",
                                  base_period = "varying",
                                  bstrap = T,
                                  clustervars = "ID",
                                  control_group = "nevertreated",
                                  #xformla = ~Population+CO2,
                                  data = did_master_cities_panel_pm
)
simple.attgt_retirement50k_var_pm <- aggte(attgt_retirement50k_var_pm, type = "simple")
dynamic.attgt_retirement50k_var_pm <- aggte(attgt_retirement50k_var_pm, type = "dynamic")
cohort.attgt_retirement50k_var_pm <- aggte(attgt_retirement50k_var_pm, type = "group")

attgt_retirement50k_var_cov_pm <- att_gt(yname = "Rates_PM",
                                      tname = "year",
                                      idname = "ID",
                                      gname = "Retired.year50k",
                                      base_period = "varying",
                                      bstrap = T,
                                      clustervars = "ID",
                                      control_group = "nevertreated",
                                      xformla = ~Population+CO2,
                                      data = did_master_cities_panel_pm
)
simple.attgt_retirement50k_var_cov_pm <- aggte(attgt_retirement50k_var_cov_pm, type = "simple")
dynamic.attgt_retirement50k_var_cov_pm <- aggte(attgt_retirement50k_var_cov_pm, type = "dynamic")
cohort.attgt_retirement50k_var_cov_pm_rate_pm <- aggte(attgt_retirement50k_var_cov_pm, type = "group")

attgt_retirement50k_uni_pm <- att_gt(yname = "Rates_PM",
                                  tname = "year",
                                  idname = "ID",
                                  gname = "Retired.year50k",
                                  base_period = "universal",
                                  bstrap = T,
                                  clustervars = "ID",
                                  control_group = "nevertreated",
                                  #xformla = ~Population+CO2,
                                  data = did_master_cities_panel_pm
)
simple.attgt_retirement50k_uni_pm <- aggte(attgt_retirement50k_uni_pm, type = "simple")
dynamic.attgt_retirement50k_uni_pm <- aggte(attgt_retirement50k_uni_pm, type = "dynamic")
cohort.attgt_retirement50k_uni_pm <- aggte(attgt_retirement50k_uni_pm, type = "group")

attgt_retirement50k_uni_cov_pm <- att_gt(yname = "Rates_PM",
                                      tname = "year",
                                      idname = "ID",
                                      gname = "Retired.year50k",
                                      base_period = "universal",
                                      bstrap = T,
                                      clustervars = "ID",
                                      control_group = "nevertreated",
                                      xformla = ~Population+CO2,
                                      data = did_master_cities_panel_pm
)
simple.attgt_retirement50k_uni_cov_pm <- aggte(attgt_retirement50k_uni_cov_pm, type = "simple")
dynamic.attgt_retirement50k_uni_cov_pm <- aggte(attgt_retirement50k_uni_cov_pm, type = "dynamic")
cohort.attgt_retirement50k_uni_cov_pm <- aggte(attgt_retirement50k_uni_cov_pm, type = "group")

mod1_pm_retirement50k_var_pm <- tidy(dynamic.attgt_retirement50k_var_pm) %>% mutate(model = "50k from Retired Coal Unit (Mod1)")
mod1_pm_retirement50k_var_cov_pm <- tidy(dynamic.attgt_retirement50k_var_cov_pm) %>% mutate(model = "50k from Retired Coal Unit (Mod2)")
mod1_pm_retirement50k_uni_pm <- tidy(dynamic.attgt_retirement50k_uni_pm) %>% mutate(model = "50k from Retired Coal Unit (Mod3)")
mod1_pm_retirement50k_uni_cov_pm <- tidy(dynamic.attgt_retirement50k_uni_cov_pm) %>% mutate(model = "50k from Retired Coal Unit (Mod4)")

#The models above test the robustness of the findings to using different ways to estimate the pre-treatment trends and the inclusion of covariates

mod1_pm_retirement <- bind_rows(mod1_pm_retirement50k_var_pm, mod1_pm_retirement50k_var_cov_pm,
                                mod1_pm_retirement50k_uni_pm, mod1_pm_retirement50k_uni_cov_pm)

pdf(file = "results/fig10_pm_main.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot(mod1_pm_retirement, aes(x = as.factor(event.time), y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange(aes(color = model, group = model), position = position_dodge(width=0.3)) +
  xlab("Time to Coal Power Unit Retirement") +
  ylab("Coefficient") +
  ggtitle("Average Effect by Length of Exposure") +
  theme_classic() + scale_color_brewer(palette = "Set1")  + 
  geom_vline(xintercept = "0", linetype = "twodash", color = "grey50", size=0.5) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", size=0.3)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure 16: Event study analysis for mortality rates associated with NO2.                                           #
#--------------------------------------------------------------------------------------------------------------------#

attgt_retirement50k_var_no2 <- att_gt(yname = "Rates_NO2",
                                  tname = "year",
                                  idname = "ID",
                                  gname = "Retired.year50k",
                                  base_period = "varying",
                                  bstrap = T,
                                  clustervars = "ID",
                                  control_group = "nevertreated",
                                  #xformla = ~Population+CO2,
                                  data = did_master_cities_panel_no2
)
simple.attgt_retirement50k_var_no2 <- aggte(attgt_retirement50k_var_no2, type = "simple")
dynamic.attgt_retirement50k_var_no2 <- aggte(attgt_retirement50k_var_no2, type = "dynamic")
cohort.attgt_retirement50k_var_no2_rate <- aggte(attgt_retirement50k_var_no2, type = "group")

attgt_retirement50k_uni_no2 <- att_gt(yname = "Rates_NO2",
                                  tname = "year",
                                  idname = "ID",
                                  gname = "Retired.year50k",
                                  base_period = "universal",
                                  bstrap = T,
                                  clustervars = "ID",
                                  control_group = "nevertreated",
                                  #xformla = ~Population+CO2,
                                  data = did_master_cities_panel_no2
)
simple.attgt_retirement50k_uni_no2 <- aggte(attgt_retirement50k_uni_no2, type = "simple")
dynamic.attgt_retirement50k_uni_no2 <- aggte(attgt_retirement50k_uni_no2, type = "dynamic")

attgt_retirement50k_var_cov_no2 <- att_gt(yname = "Rates_NO2",
                                      tname = "year",
                                      idname = "ID",
                                      gname = "Retired.year50k",
                                      base_period = "varying",
                                      bstrap = T,
                                      clustervars = "ID",
                                      control_group = "nevertreated",
                                      xformla = ~Population+CO2,
                                      data = did_master_cities_panel_no2
)
simple.attgt_retirement50k_var_cov_no2 <- aggte(attgt_retirement50k_var_cov_no2, type = "simple")
dynamic.attgt_retirement50k_var_cov_no2 <- aggte(attgt_retirement50k_var_cov_no2, type = "dynamic")

attgt_retirement50k_uni_cov_no2 <- att_gt(yname = "Rates_NO2",
                                      tname = "year",
                                      idname = "ID",
                                      gname = "Retired.year50k",
                                      base_period = "universal",
                                      bstrap = T,
                                      clustervars = "ID",
                                      control_group = "nevertreated",
                                      xformla = ~Population+CO2,
                                      data = did_master_cities_panel_no2
)
simple.attgt_retirement50k_uni_cov_no2 <- aggte(attgt_retirement50k_uni_cov_no2, type = "simple")
dynamic.attgt_retirement50k_uni_cov_no2 <- aggte(attgt_retirement50k_uni_cov_no2, type = "dynamic")

mod1_no2_retirement50k_var <- tidy(dynamic.attgt_retirement50k_var_no2) %>% mutate(model = "50k from Retired Coal Unit (Mod1)")
mod1_no2_retirement50k_var_cov <- tidy(dynamic.attgt_retirement50k_var_cov_no2) %>% mutate(model = "50k from Retired Coal Unit (Mod2)")
mod1_no2_retirement50k_uni <- tidy(dynamic.attgt_retirement50k_uni_no2) %>% mutate(model = "50k from Retired Coal Unit (Mod3)")
mod1_no2_retirement50k_uni_cov <- tidy(dynamic.attgt_retirement50k_uni_cov_no2) %>% mutate(model = "50k from Retired Coal Unit (Mod4)")

mod1_no2_retirement <- bind_rows(mod1_no2_retirement50k_var, mod1_no2_retirement50k_var_cov,
                                mod1_no2_retirement50k_uni, mod1_no2_retirement50k_uni_cov)

pdf(file = "results/fig11_no2_main.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot(mod1_no2_retirement, aes(x = as.factor(event.time), y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange(aes(color = model, group = model), position = position_dodge(width=0.3)) +
  xlab("Time to Coal Power Unit Retirement") +
  ylab("Coefficient") +
  ggtitle("Average Effect by Length of Exposure") +
  theme_classic() + scale_color_brewer(palette = "Set1")  + 
  geom_vline(xintercept = "0", linetype = "twodash", color = "grey50", size=0.5) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", size=0.3)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure 17: Heterogeneous effects of phasing out coal power plants on public health outcomes (Mortality rates-PM)   #                                      #
#--------------------------------------------------------------------------------------------------------------------#

#Left-Panel

attgt_retirement100k_var_pm <- att_gt(yname = "Rates_PM",
                                   tname = "year",
                                   idname = "ID",
                                   gname = "Retired.year100k",
                                   base_period = "varying",
                                   bstrap = T,
                                   clustervars = "ID",
                                   control_group = "nevertreated",
                                   #xformla = ~Population+CO2,
                                   data = did_master_cities_panel_pm
)
simple.attgt_retirement100k_var_pm <- aggte(attgt_retirement100k_var_pm, type = "simple")
dynamic.attgt_retirement100k_var_pm <- aggte(attgt_retirement100k_var_pm, type = "dynamic")
cohort.attgt_retirement100k_var_pm <- aggte(attgt_retirement100k_var_pm, type = "group")

attgt_retirement150k_var_pm <- att_gt(yname = "Rates_PM",
                                   tname = "year",
                                   idname = "ID",
                                   gname = "Retired.year150k",
                                   base_period = "varying",
                                   bstrap = T,
                                   clustervars = "ID",
                                   control_group = "nevertreated",
                                   #xformla = ~Population+CO2,
                                   data = did_master_cities_panel_pm
)
simple.attgt_retirement150k_var_pm <- aggte(attgt_retirement150k_var_pm, type = "simple")
dynamic.attgt_retirement150k_var_pm <- aggte(attgt_retirement150k_var_pm, type = "dynamic")

mod1_pm_retirement50k_var <- tidy(dynamic.attgt_retirement50k_var_pm) %>% mutate(model = "50k from Retired Coal Unit (Mod1)")
mod1_pm_retirement100k_var <- tidy(dynamic.attgt_retirement100k_var_pm) %>% mutate(model = "100k from Retired Coal Unit")
mod1_pm_retirement150k_var <- tidy(dynamic.attgt_retirement150k_var_pm) %>% mutate(model = "150k from Retired Coal Unit")

mod1_pm_retirement_robustness1_dist <- bind_rows(mod1_pm_retirement50k_var, mod1_pm_retirement100k_var, mod1_pm_retirement150k_var)

pdf(file = "results/fig10b_pm_distance.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot(mod1_pm_retirement_robustness1_dist, aes(x = as.factor(event.time), y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange(aes(color = model, group = model), position = position_dodge(width=0.3)) +
  xlab("Time to Coal Power Unit Retirement") +
  ylab("Coefficient") +
  ggtitle("Average Effect by Length of Exposure") +
  theme_classic() + scale_color_brewer(palette = "Set1")  + 
  geom_vline(xintercept = "0", linetype = "twodash", color = "grey50", size=0.5) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", size=0.3)
dev.off()

#Right-side panel

attgt_retirement50k_high_var_pm <- att_gt(yname = "Rates_PM",
                                       tname = "year",
                                       idname = "ID",
                                       gname = "Retired.year50k_high",
                                       base_period = "varying",
                                       bstrap = T,
                                       clustervars = "ID",
                                       control_group = "nevertreated",
                                       #xformla = ~Population+CO2,
                                       data = did_master_cities_panel_pm
)
simple.attgt_retirement50k_high_var_pm <- aggte(attgt_retirement50k_high_var_pm, type = "simple")
dynamic.attgt_retirement50k_high_var_pm <- aggte(attgt_retirement50k_high_var_pm, type = "dynamic")
cohort.attgt_retirement50k_high_var_pm_rate <- aggte(attgt_retirement50k_high_var_pm, type = "group")

attgt_retirement50k_low_var_pm <- att_gt(yname = "Rates_PM",
                                      tname = "year",
                                      idname = "ID",
                                      gname = "Retired.year50k_low",
                                      base_period = "varying",
                                      bstrap = T,
                                      clustervars = "ID",
                                      control_group = "nevertreated",
                                      #xformla = ~Population+CO2,
                                      data = did_master_cities_panel_pm
)
simple.attgt_retirement50k_low_var_pm <- aggte(attgt_retirement50k_low_var_pm, type = "simple")
dynamic.attgt_retirement50k_low_var_pm <- aggte(attgt_retirement50k_low_var_pm, type = "dynamic")
cohort.attgt_retirement50k_low_var_pm_rate <- aggte(attgt_retirement50k_low_var_pm, type = "group")

mod1_pm_retirement50k_high_var <- tidy(dynamic.attgt_retirement50k_high_var_pm) %>% mutate(model = "50k from Retired Coal Unit (High Capacity MW)")
mod1_pm_retirement50k_low_var <- tidy(dynamic.attgt_retirement50k_low_var_pm) %>% mutate(model = "50k from Retired Coal Unit (Low Capacity MW)")

mod1_pm_retirement_robustness2_capacity <- bind_rows(mod1_pm_retirement50k_high_var, mod1_pm_retirement50k_low_var)

pdf(file = "results/fig10c_pm_capacity.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot(mod1_pm_retirement_robustness2_capacity, aes(x = as.factor(event.time), y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange(aes(color = model, group = model), position = position_dodge(width=0.3)) +
  xlab("Time to Coal Power Unit Retirement") +
  ylab("Coefficient") +
  ggtitle("Average Effect by Length of Exposure") +
  theme_classic() + scale_color_brewer(palette = "Set1")  + 
  geom_vline(xintercept = "0", linetype = "twodash", color = "grey50", size=0.5) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", size=0.3)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure 18: Heterogeneous effects of phasing out coal power plants on public health outcomes (Mortality rates-NO2)  #                                          #
#--------------------------------------------------------------------------------------------------------------------#

#Left-Panel

attgt_retirement100k_var_no2 <- att_gt(yname = "Rates_NO2",
                                   tname = "year",
                                   idname = "ID",
                                   gname = "Retired.year100k",
                                   base_period = "varying",
                                   bstrap = T,
                                   clustervars = "ID",
                                   control_group = "nevertreated",
                                   #xformla = ~Population+CO2,
                                   data = did_master_cities_panel_no2
)
simple.attgt_retirement100k_var_no2 <- aggte(attgt_retirement100k_var_no2, type = "simple")
dynamic.attgt_retirement100k_var_no2 <- aggte(attgt_retirement100k_var_no2, type = "dynamic")

attgt_retirement150k_var_no2 <- att_gt(yname = "Rates_NO2",
                                   tname = "year",
                                   idname = "ID",
                                   gname = "Retired.year150k",
                                   base_period = "varying",
                                   bstrap = T,
                                   clustervars = "ID",
                                   control_group = "nevertreated",
                                   #xformla = ~Population+CO2,
                                   data = did_master_cities_panel_no2
)
simple.attgt_retirement150k_var_no2 <- aggte(attgt_retirement150k_var_no2, type = "simple")
dynamic.attgt_retirement150k_var_no2 <- aggte(attgt_retirement150k_var_no2, type = "dynamic")

mod1_no2_retirement100k_var <- tidy(dynamic.attgt_retirement100k_var_no2) %>% mutate(model = "100k from Retired Coal Unit")
mod1_no2_retirement150k_var <- tidy(dynamic.attgt_retirement150k_var_no2) %>% mutate(model = "150k from Retired Coal Unit")

mod1_no2_retirement_robustness1_dist <- bind_rows(mod1_no2_retirement50k_var, mod1_no2_retirement100k_var, mod1_no2_retirement150k_var)

pdf(file = "results/fig11b_no2_distance.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot(mod1_no2_retirement_robustness1_dist, aes(x = as.factor(event.time), y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange(aes(color = model, group = model), position = position_dodge(width=0.3)) +
  xlab("Time to Coal Power Unit Retirement") +
  ylab("Coefficient") +
  ggtitle("Average Effect by Length of Exposure") +
  theme_classic() + scale_color_brewer(palette = "Set1")  + 
  geom_vline(xintercept = "0", linetype = "twodash", color = "grey50", size=0.5) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", size=0.3)
dev.off()

#Right Panel

attgt_retirement50k_high_var_no2 <- att_gt(yname = "Rates_NO2",
                                       tname = "year",
                                       idname = "ID",
                                       gname = "Retired.year50k_high",
                                       base_period = "varying",
                                       bstrap = T,
                                       clustervars = "ID",
                                       control_group = "nevertreated",
                                       #xformla = ~Population+CO2,
                                       data = did_master_cities_panel_no2)
simple.attgt_retirement50k_high_var_no2 <- aggte(attgt_retirement50k_high_var_no2, type = "simple")
dynamic.attgt_retirement50k_high_var_no2 <- aggte(attgt_retirement50k_high_var_no2, type = "dynamic")
cohort.attgt_retirement50k_high_var_no2_rate <- aggte(attgt_retirement50k_high_var_no2, type = "group")

attgt_retirement50k_low_var_no2 <- att_gt(yname = "Rates_NO2",
                                      tname = "year",
                                      idname = "ID",
                                      gname = "Retired.year50k_low",
                                      base_period = "varying",
                                      bstrap = T,
                                      clustervars = "ID",
                                      control_group = "nevertreated",
                                      #xformla = ~Population+CO2,
                                      data = did_master_cities_panel_no2)
simple.attgt_retirement50k_low_var_no2 <- aggte(attgt_retirement50k_low_var_no2, type = "simple")
dynamic.attgt_retirement50k_low_var_no2 <- aggte(attgt_retirement50k_low_var_no2, type = "dynamic")
cohort.attgt_retirement50k_low_var_no2_rate <- aggte(attgt_retirement50k_low_var_no2, type = "group")

mod1_no2_retirement50k_high_var <- tidy(dynamic.attgt_retirement50k_high_var_no2) %>% mutate(model = "50k from Retired Coal Unit (High Capacity MW)")
mod1_no2_retirement50k_low_var <- tidy(dynamic.attgt_retirement50k_low_var_no2) %>% mutate(model = "50k from Retired Coal Unit (Low Capacity MW)")

mod1_no2_retirement_robustness2_capacity <- bind_rows(mod1_no2_retirement50k_high_var, mod1_no2_retirement50k_low_var)

pdf(file = "results/fig11c_no2_capacity.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot(mod1_no2_retirement_robustness2_capacity, aes(x = as.factor(event.time), y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange(aes(color = model, group = model), position = position_dodge(width=0.3)) +
  xlab("Time to Coal Power Unit Retirement") +
  ylab("Coefficient") +
  ggtitle("Average Effect by Length of Exposure") +
  theme_classic() + scale_color_brewer(palette = "Set1")  + 
  geom_vline(xintercept = "0", linetype = "twodash", color = "grey50", size=0.5) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", size=0.3)
dev.off()

##############################################
####### FIGURES SUPPORTING INFORMATION #######
##############################################

#Creating data for PM2.5

did_master_cities_panel_pm_subset <- master_cities_panel %>%
  filter(close_retired50k == 1 | close_operating100k == 1) %>%
  mutate(Retired.year50k = ifelse(is.na(Retired.year50k), 0, Retired.year50k),
         Retired.year100k = ifelse(is.na(Retired.year100k), 0, Retired.year100k),
         Retired.year150k = ifelse(is.na(Retired.year150k), 0, Retired.year150k),
         Retired.year450k = ifelse(is.na(Retired.year450k), 0, Retired.year450k),
         Retired.year50k_high = ifelse(is.na(Retired.year50k_high), 0, Retired.year50k_high),
         Retired.year50k_low = ifelse(is.na(Retired.year50k_low), 0, Retired.year50k_low)) %>%
  dplyr::select(ID, year, Rates_PM, Retired.year50k, Retired.year100k, 
                Retired.year150k, Retired.year450k, Population, CO2, political_regime,
                Retired.year50k_high, Retired.year50k_low) %>% na.omit()

#Creating data for NO2

did_master_cities_panel_no2_subset <- master_cities_panel %>%
  filter(close_retired50k == 1 | close_operating100k == 1) %>%
  mutate(Retired.year50k = ifelse(is.na(Retired.year50k), 0, Retired.year50k),
         Retired.year100k = ifelse(is.na(Retired.year100k), 0, Retired.year100k),
         Retired.year150k = ifelse(is.na(Retired.year150k), 0, Retired.year150k),
         Retired.year450k = ifelse(is.na(Retired.year450k), 0, Retired.year450k),
         Retired.year50k_high = ifelse(is.na(Retired.year50k_high), 0, Retired.year50k_high),
         Retired.year50k_low = ifelse(is.na(Retired.year50k_low), 0, Retired.year50k_low)) %>%
  dplyr::select(ID, year, Rates_NO2, Retired.year50k, Retired.year100k, 
                Retired.year150k, Retired.year450k, Population, CO2, political_regime,
                Retired.year50k_high, Retired.year50k_low) %>% na.omit()

#Creating data for PM

did_master_cities_panel_pm_emissions_subset <- master_cities_panel %>%
  filter(close_retired50k == 1 | close_operating100k == 1) %>%
  mutate(Retired.year50k = ifelse(is.na(Retired.year50k), 0, Retired.year50k),
         Retired.year100k = ifelse(is.na(Retired.year100k), 0, Retired.year100k),
         Retired.year150k = ifelse(is.na(Retired.year150k), 0, Retired.year150k),
         Retired.year450k = ifelse(is.na(Retired.year450k), 0, Retired.year450k),
         Retired.year50k_high = ifelse(is.na(Retired.year50k_high), 0, Retired.year50k_high),
         Retired.year50k_low = ifelse(is.na(Retired.year50k_low), 0, Retired.year50k_low)) %>%
  dplyr::select(ID, year, PM, Retired.year50k, Retired.year100k, 
                Retired.year150k, Retired.year450k, Population, CO2, political_regime,
                Retired.year50k_high, Retired.year50k_low) %>% na.omit()

did_master_cities_panel_no2_emissions_subset <- master_cities_panel %>%
  filter(close_retired50k == 1 | close_operating100k == 1) %>%
  mutate(Retired.year50k = ifelse(is.na(Retired.year50k), 0, Retired.year50k),
         Retired.year100k = ifelse(is.na(Retired.year100k), 0, Retired.year100k),
         Retired.year150k = ifelse(is.na(Retired.year150k), 0, Retired.year150k),
         Retired.year450k = ifelse(is.na(Retired.year450k), 0, Retired.year450k),
         Retired.year50k_high = ifelse(is.na(Retired.year50k_high), 0, Retired.year50k_high),
         Retired.year50k_low = ifelse(is.na(Retired.year50k_low), 0, Retired.year50k_low)) %>%
  dplyr::select(ID, year, NO2, Retired.year50k, Retired.year100k, 
                Retired.year150k, Retired.year450k, Population, CO2, political_regime,
                Retired.year50k_high, Retired.year50k_low) %>% na.omit()

did_master_cities_panel_co2_emissions_subset <- master_cities_panel %>%
  filter(close_retired50k == 1 | close_operating100k == 1) %>%
  mutate(Retired.year50k = ifelse(is.na(Retired.year50k), 0, Retired.year50k),
         Retired.year100k = ifelse(is.na(Retired.year100k), 0, Retired.year100k),
         Retired.year150k = ifelse(is.na(Retired.year150k), 0, Retired.year150k),
         Retired.year450k = ifelse(is.na(Retired.year450k), 0, Retired.year450k),
         Retired.year50k_high = ifelse(is.na(Retired.year50k_high), 0, Retired.year50k_high),
         Retired.year50k_low = ifelse(is.na(Retired.year50k_low), 0, Retired.year50k_low)) %>%
  dplyr::select(ID, year, NO2, Retired.year50k, Retired.year100k, 
                Retired.year150k, Retired.year450k, Population, CO2, political_regime,
                Retired.year50k_high, Retired.year50k_low) %>% na.omit() %>% filter(CO2 > 0) %>%
  mutate(logco2 = log(CO2))

#--------------------------------------------------------------------------------------------------------------------#
# Figure S1: Staggered difference-in-differences for PM2.5 Rates                                                     # 
#--------------------------------------------------------------------------------------------------------------------#

attgt_retirement50k_var_subset <- att_gt(yname = "Rates_PM",
                                  tname = "year",
                                  idname = "ID",
                                  gname = "Retired.year50k",
                                  base_period = "varying",
                                  bstrap = T,
                                  clustervars = "ID",
                                  control_group = "nevertreated",
                                  #xformla = ~Population+CO2,
                                  data = did_master_cities_panel_pm_subset
)
simple.attgt_retirement50k_var_subset <- aggte(attgt_retirement50k_var_subset, type = "simple")
dynamic.attgt_retirement50k_var_subset <- aggte(attgt_retirement50k_var_subset, type = "dynamic")

attgt_retirement50k_high_var_subset <- att_gt(yname = "Rates_PM",
                                       tname = "year",
                                       idname = "ID",
                                       gname = "Retired.year50k_high",
                                       base_period = "varying",
                                       bstrap = T,
                                       clustervars = "ID",
                                       control_group = "nevertreated",
                                       #xformla = ~Population+CO2,
                                       data = did_master_cities_panel_pm_subset
)
simple.attgt_retirement50k_high_var_subset <- aggte(attgt_retirement50k_high_var_subset, type = "simple")
dynamic.attgt_retirement50k_high_var_subset <- aggte(attgt_retirement50k_high_var_subset, type = "dynamic")

attgt_retirement50k_low_var_subset <- att_gt(yname = "Rates_PM",
                                      tname = "year",
                                      idname = "ID",
                                      gname = "Retired.year50k_low",
                                      base_period = "varying",
                                      bstrap = T,
                                      clustervars = "ID",
                                      control_group = "nevertreated",
                                      #xformla = ~Population+CO2,
                                      data = did_master_cities_panel_pm_subset
)
simple.attgt_retirement50k_low_var_subset <- aggte(attgt_retirement50k_low_var_subset, type = "simple")
dynamic.attgt_retirement50k_low_var_subset <- aggte(attgt_retirement50k_low_var_subset, type = "dynamic")

mod1_pm_retirement50k_var_subset <- tidy(dynamic.attgt_retirement50k_var_subset) %>% mutate(model = "50k from Retired Coal Unit (All)")
mod1_pm_retirement50k_high_var_subset <- tidy(dynamic.attgt_retirement50k_high_var_subset) %>% mutate(model = "50k from Retired Coal Unit (High Capacity MW)")
mod1_pm_retirement50k_low_var_subset <- tidy(dynamic.attgt_retirement50k_low_var_subset) %>% mutate(model = "50k from Retired Coal Unit (Low Capacity MW)")

mod1_pm_retirement_subset <- bind_rows(mod1_pm_retirement50k_var_subset, mod1_pm_retirement50k_high_var_subset, mod1_pm_retirement50k_low_var_subset)

pdf(file = "si_figures/figSI1_pm_main.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot(mod1_pm_retirement_subset, aes(x = as.factor(event.time), y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange(aes(color = model, group = model), position = position_dodge(width=0.3)) +
  xlab("Time to Coal Power Unit Retirement") +
  ylab("Coefficient") +
  ggtitle("Average Effect by Length of Exposure") +
  theme_classic() + scale_color_brewer(palette = "Set1")  + 
  geom_vline(xintercept = "0", linetype = "twodash", color = "grey50", size=0.5) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", size=0.3)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S2: Staggered difference-in-differences for NO2 Rates.                                                      # 
#--------------------------------------------------------------------------------------------------------------------#

attgt_retirement50k_var_subset <- att_gt(yname = "Rates_NO2",
                                  tname = "year",
                                  idname = "ID",
                                  gname = "Retired.year50k",
                                  base_period = "varying",
                                  bstrap = T,
                                  clustervars = "ID",
                                  control_group = "nevertreated",
                                  #xformla = ~Population+CO2,
                                  data = did_master_cities_panel_no2_subset
)
simple.attgt_retirement50k_var_subset <- aggte(attgt_retirement50k_var_subset, type = "simple")
dynamic.attgt_retirement50k_var_subset <- aggte(attgt_retirement50k_var_subset, type = "dynamic")

attgt_retirement50k_high_var_subset <- att_gt(yname = "Rates_NO2",
                                       tname = "year",
                                       idname = "ID",
                                       gname = "Retired.year50k_high",
                                       base_period = "varying",
                                       bstrap = T,
                                       clustervars = "ID",
                                       control_group = "nevertreated",
                                       #xformla = ~Population+CO2,
                                       data = did_master_cities_panel_no2_subset
)
simple.attgt_retirement50k_high_var_subset <- aggte(attgt_retirement50k_high_var_subset, type = "simple")
dynamic.attgt_retirement50k_high_var_subset <- aggte(attgt_retirement50k_high_var_subset, type = "dynamic")

attgt_retirement50k_low_var_subset <- att_gt(yname = "Rates_NO2",
                                      tname = "year",
                                      idname = "ID",
                                      gname = "Retired.year50k_low",
                                      base_period = "varying",
                                      bstrap = T,
                                      clustervars = "ID",
                                      control_group = "nevertreated",
                                      #xformla = ~Population+CO2,
                                      data = did_master_cities_panel_no2_subset
)
simple.attgt_retirement50k_low_var_subset <- aggte(attgt_retirement50k_low_var_subset, type = "simple")
dynamic.attgt_retirement50k_low_var_subset <- aggte(attgt_retirement50k_low_var_subset, type = "dynamic")

mod1_no2_retirement50k_var_subset <- tidy(dynamic.attgt_retirement50k_var_subset) %>% mutate(model = "50k from Retired Coal Unit (All)")
mod1_no2_retirement50k_high_var_subset <- tidy(dynamic.attgt_retirement50k_high_var_subset) %>% mutate(model = "50k from Retired Coal Unit (High Capacity MW)")
mod1_no2_retirement50k_low_var_subset <- tidy(dynamic.attgt_retirement50k_low_var_subset) %>% mutate(model = "50k from Retired Coal Unit (Low Capacity MW)")

mod1_no2_retirement_subset <- bind_rows(mod1_no2_retirement50k_var_subset, mod1_no2_retirement50k_high_var_subset, mod1_no2_retirement50k_low_var_subset)

pdf(file = "si_figures/figSI2_no2_main.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot(mod1_no2_retirement_subset, aes(x = as.factor(event.time), y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange(aes(color = model, group = model), position = position_dodge(width=0.3)) +
  xlab("Time to Coal Power Unit Retirement") +
  ylab("Coefficient") +
  ggtitle("Average Effect by Length of Exposure") +
  theme_classic() + scale_color_brewer(palette = "Set1")  + 
  geom_vline(xintercept = "0", linetype = "twodash", color = "grey50", size=0.5) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", size=0.3)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S3: Staggered differences-in-differences, cohort aggregation (PM2.5 Rates).                                 # 
#--------------------------------------------------------------------------------------------------------------------#
pdf(file = "si_figures/figSI15cohort_pm_rates.main.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggdid(cohort.attgt_retirement50k_var_pm)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S4: Staggered differences-in-differences, cohort aggregation (PM2.5 Rates).                                 # 
#--------------------------------------------------------------------------------------------------------------------#

pdf(file = "si_figures/figSI16cohort_pm_rates.high.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggdid(cohort.attgt_retirement50k_high_var_pm_rate)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S5: Staggered differences-in-differences, cohort aggregation (PM2.5 Rates).                                 # 
#--------------------------------------------------------------------------------------------------------------------#

pdf(file = "si_figures/figSI17cohort_pm_rates.low.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggdid(cohort.attgt_retirement50k_low_var_pm_rate)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S6: Staggered differences-in-differences, cohort aggregation (NO2 Rates).                                   # 
#--------------------------------------------------------------------------------------------------------------------#

pdf(file = "si_figures/figSI18cohort_no2_rates.main.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggdid(cohort.attgt_retirement50k_var_no2_rate)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S7: Staggered differences-in-differences, cohort aggregation (NO2 Rates).                                   # 
#--------------------------------------------------------------------------------------------------------------------#

pdf(file = "si_figures/figSI19cohort_no2_rates.high.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggdid(cohort.attgt_retirement50k_high_var_no2_rate)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S8: Staggered differences-in-differences, cohort aggregation (NO2 Rates).                                   # 
#--------------------------------------------------------------------------------------------------------------------#

pdf(file = "si_figures/figSI20cohort_no2_rates.low.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggdid(cohort.attgt_retirement50k_low_var_no2_rate)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S9: Staggered differences-in-differences, event study (PM2.5 Rates).                                        # 
#--------------------------------------------------------------------------------------------------------------------#

#Creating data for PM2.5 Rates

did_master_cities_panel_pm_robustness <- master_cities_panel_robustness %>%
  mutate(Retired.year50k = ifelse(is.na(Retired.year50k), 0, Retired.year50k),
         Retired.year100k = ifelse(is.na(Retired.year100k), 0, Retired.year100k),
         Retired.year150k = ifelse(is.na(Retired.year150k), 0, Retired.year150k),
         Retired.year450k = ifelse(is.na(Retired.year450k), 0, Retired.year450k),
         Retired.year50k_high = ifelse(is.na(Retired.year50k_high), 0, Retired.year50k_high),
         Retired.year50k_low = ifelse(is.na(Retired.year50k_low), 0, Retired.year50k_low)) %>%
  dplyr::select(ID, year, Rates_PM, Retired.year50k, Retired.year100k, 
                Retired.year150k, Retired.year450k, Population, CO2, political_regime,
                Retired.year50k_high, Retired.year50k_low) %>% na.omit()

#Creating data for NO2 Rates

did_master_cities_panel_no2_robustness <- master_cities_panel_robustness %>%
  mutate(Retired.year50k = ifelse(is.na(Retired.year50k), 0, Retired.year50k),
         Retired.year100k = ifelse(is.na(Retired.year100k), 0, Retired.year100k),
         Retired.year150k = ifelse(is.na(Retired.year150k), 0, Retired.year150k),
         Retired.year450k = ifelse(is.na(Retired.year450k), 0, Retired.year450k),
         Retired.year50k_high = ifelse(is.na(Retired.year50k_high), 0, Retired.year50k_high),
         Retired.year50k_low = ifelse(is.na(Retired.year50k_low), 0, Retired.year50k_low)) %>%
  dplyr::select(ID, year, Rates_NO2, Retired.year50k, Retired.year100k, 
                Retired.year150k, Retired.year450k, Population, CO2, political_regime,
                Retired.year50k_high, Retired.year50k_low) %>% na.omit()

attgt_retirement50k_var_retonly <- att_gt(yname = "Rates_PM",
                                  tname = "year",
                                  idname = "ID",
                                  gname = "Retired.year50k",
                                  base_period = "varying",
                                  bstrap = T,
                                  clustervars = "ID",
                                  control_group = "nevertreated",
                                  #xformla = ~Population+CO2,
                                  data = did_master_cities_panel_pm_robustness
)
simple.attgt_retirement50k_var_retonly <- aggte(attgt_retirement50k_var_retonly, type = "simple")
dynamic.attgt_retirement50k_var_retonly <- aggte(attgt_retirement50k_var_retonly, type = "dynamic")
cohort.attgt_retirement50k_var_retonly <- aggte(attgt_retirement50k_var_retonly, type = "group")

attgt_retirement50k_uni_retonly <- att_gt(yname = "Rates_PM",
                                  tname = "year",
                                  idname = "ID",
                                  gname = "Retired.year50k",
                                  base_period = "universal",
                                  bstrap = T,
                                  clustervars = "ID",
                                  control_group = "nevertreated",
                                  #xformla = ~Population+CO2,
                                  data = did_master_cities_panel_pm_robustness
)
simple.attgt_retirement50k_uni_retonly <- aggte(attgt_retirement50k_uni_retonly, type = "simple")
dynamic.attgt_retirement50k_uni_retonly <- aggte(attgt_retirement50k_uni_retonly, type = "dynamic")
cohort.attgt_retirement50k_uni_retonly <- aggte(attgt_retirement50k_uni_retonly, type = "group")

attgt_retirement50k_var_cov_retonly <- att_gt(yname = "Rates_PM",
                                      tname = "year",
                                      idname = "ID",
                                      gname = "Retired.year50k",
                                      base_period = "varying",
                                      bstrap = T,
                                      clustervars = "ID",
                                      control_group = "nevertreated",
                                      xformla = ~Population+CO2,
                                      data = did_master_cities_panel_pm_robustness
)
simple.attgt_retirement50k_var_cov_retonly <- aggte(attgt_retirement50k_var_cov_retonly, type = "simple")
dynamic.attgt_retirement50k_var_cov_retonly <- aggte(attgt_retirement50k_var_cov_retonly, type = "dynamic")
cohort.attgt_retirement50k_var_cov_pm_rate_retonly <- aggte(attgt_retirement50k_var_cov_retonly, type = "group")

attgt_retirement50k_uni_cov_retonly <- att_gt(yname = "Rates_PM",
                                      tname = "year",
                                      idname = "ID",
                                      gname = "Retired.year50k",
                                      base_period = "universal",
                                      bstrap = T,
                                      clustervars = "ID",
                                      control_group = "nevertreated",
                                      xformla = ~Population+CO2,
                                      data = did_master_cities_panel_pm_robustness
)
simple.attgt_retirement50k_uni_cov_retonly <- aggte(attgt_retirement50k_uni_cov_retonly, type = "simple")
dynamic.attgt_retirement50k_uni_cov_retonly <- aggte(attgt_retirement50k_uni_cov_retonly, type = "dynamic")
cohort.attgt_retirement50k_uni_cov_retonly <- aggte(attgt_retirement50k_uni_cov_retonly, type = "group")

mod1_pm_retirement50k_var_retonly <- tidy(dynamic.attgt_retirement50k_var_retonly) %>% mutate(model = "50k from Retired Coal Unit (Mod1)")
mod1_pm_retirement50k_var_cov_retonly <- tidy(dynamic.attgt_retirement50k_var_cov_retonly) %>% mutate(model = "50k from Retired Coal Unit (Mod2)")
mod1_pm_retirement50k_uni_retonly <- tidy(dynamic.attgt_retirement50k_uni_retonly) %>% mutate(model = "50k from Retired Coal Unit (Mod3)")
mod1_pm_retirement50k_uni_cov_retonly <- tidy(dynamic.attgt_retirement50k_uni_cov_retonly) %>% mutate(model = "50k from Retired Coal Unit (Mod4)")

mod1_pm_retirement_retonly <- bind_rows(mod1_pm_retirement50k_var_retonly, mod1_pm_retirement50k_var_cov_retonly,
                                mod1_pm_retirement50k_uni_retonly, mod1_pm_retirement50k_uni_cov_retonly)

pdf(file = "si_figures/figSI21_pm_main_onlyretired.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot(mod1_pm_retirement_retonly, aes(x = as.factor(event.time), y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange(aes(color = model, group = model), position = position_dodge(width=0.3)) +
  xlab("Time to Coal Power Unit Retirement") +
  ylab("Coefficient") +
  ggtitle("Average Effect by Length of Exposure") +
  theme_classic() + scale_color_brewer(palette = "Set1")  + 
  geom_vline(xintercept = "0", linetype = "twodash", color = "grey50", size=0.5) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", size=0.3)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S10: Staggered differences-in-differences, event study (PM2.5 Rates) by capacity.                           # 
#--------------------------------------------------------------------------------------------------------------------#

attgt_retirement50k_high_var_robustness <- att_gt(yname = "Rates_PM",
                                       tname = "year",
                                       idname = "ID",
                                       gname = "Retired.year50k_high",
                                       base_period = "varying",
                                       bstrap = T,
                                       clustervars = "ID",
                                       control_group = "nevertreated",
                                       #xformla = ~Population+CO2,
                                       data = did_master_cities_panel_pm_robustness
)
simple.attgt_retirement50k_high_var_robustness <- aggte(attgt_retirement50k_high_var_robustness, type = "simple")
dynamic.attgt_retirement50k_high_var_robustness <- aggte(attgt_retirement50k_high_var_robustness, type = "dynamic")
cohort.attgt_retirement50k_high_var_pm_rate_robustness <- aggte(attgt_retirement50k_high_var_robustness, type = "group")

attgt_retirement50k_low_var_robustness <- att_gt(yname = "Rates_PM",
                                      tname = "year",
                                      idname = "ID",
                                      gname = "Retired.year50k_low",
                                      base_period = "varying",
                                      bstrap = T,
                                      clustervars = "ID",
                                      control_group = "nevertreated",
                                      #xformla = ~Population+CO2,
                                      data = did_master_cities_panel_pm_robustness
)
simple.attgt_retirement50k_low_var_robustness <- aggte(attgt_retirement50k_low_var_robustness, type = "simple")
dynamic.attgt_retirement50k_low_var_robustness <- aggte(attgt_retirement50k_low_var_robustness, type = "dynamic")
cohort.attgt_retirement50k_low_var_pm_rate_robustness <- aggte(attgt_retirement50k_low_var_robustness, type = "group")

mod1_pm_retirement50k_high_var_robustness <- tidy(dynamic.attgt_retirement50k_high_var_robustness) %>% mutate(model = "50k from Retired Coal Unit (High Capacity MW)")
mod1_pm_retirement50k_low_var_robustness <- tidy(dynamic.attgt_retirement50k_low_var_robustness) %>% mutate(model = "50k from Retired Coal Unit (Low Capacity MW)")

mod1_pm_retirement_capacity_robustness <- bind_rows(mod1_pm_retirement50k_high_var_robustness, mod1_pm_retirement50k_low_var_robustness)

pdf(file = "si_figures/figSI22_pm_capacity_onlyretired.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot(mod1_pm_retirement_capacity_robustness, aes(x = as.factor(event.time), y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange(aes(color = model, group = model), position = position_dodge(width=0.3)) +
  xlab("Time to Coal Power Unit Retirement") +
  ylab("Coefficient") +
  ggtitle("Average Effect by Length of Exposure") +
  theme_classic() + scale_color_brewer(palette = "Set1")  + 
  geom_vline(xintercept = "0", linetype = "twodash", color = "grey50", size=0.5) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", size=0.3)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S11: Staggered differences-in-differences, event study (NO2 Rates).                                         # 
#--------------------------------------------------------------------------------------------------------------------#

attgt_retirement50k_var_robustness <- att_gt(yname = "Rates_NO2",
                                  tname = "year",
                                  idname = "ID",
                                  gname = "Retired.year50k",
                                  base_period = "varying",
                                  bstrap = T,
                                  clustervars = "ID",
                                  control_group = "nevertreated",
                                  #xformla = ~Population+CO2,
                                  data = did_master_cities_panel_no2_robustness
)
simple.attgt_retirement50k_var_robustness <- aggte(attgt_retirement50k_var_robustness, type = "simple")
dynamic.attgt_retirement50k_var_robustness <- aggte(attgt_retirement50k_var_robustness, type = "dynamic")
cohort.attgt_retirement50k_var_no2_rate_robustness <- aggte(attgt_retirement50k_var_robustness, type = "group")

attgt_retirement50k_uni_robustness <- att_gt(yname = "Rates_NO2",
                                  tname = "year",
                                  idname = "ID",
                                  gname = "Retired.year50k",
                                  base_period = "universal",
                                  bstrap = T,
                                  clustervars = "ID",
                                  control_group = "nevertreated",
                                  #xformla = ~Population+CO2,
                                  data = did_master_cities_panel_no2_robustness
)
simple.attgt_retirement50k_uni_robustness <- aggte(attgt_retirement50k_uni_robustness, type = "simple")
dynamic.attgt_retirement50k_uni_robustness <- aggte(attgt_retirement50k_uni_robustness, type = "dynamic")

attgt_retirement50k_var_cov_robustness <- att_gt(yname = "Rates_NO2",
                                      tname = "year",
                                      idname = "ID",
                                      gname = "Retired.year50k",
                                      base_period = "varying",
                                      bstrap = T,
                                      clustervars = "ID",
                                      control_group = "nevertreated",
                                      xformla = ~Population+CO2,
                                      data = did_master_cities_panel_no2_robustness
)
simple.attgt_retirement50k_var_cov_robustness <- aggte(attgt_retirement50k_var_cov_robustness, type = "simple")
dynamic.attgt_retirement50k_var_cov_robustness <- aggte(attgt_retirement50k_var_cov_robustness, type = "dynamic")

attgt_retirement50k_uni_cov_robustness <- att_gt(yname = "Rates_NO2",
                                      tname = "year",
                                      idname = "ID",
                                      gname = "Retired.year50k",
                                      base_period = "universal",
                                      bstrap = T,
                                      clustervars = "ID",
                                      control_group = "nevertreated",
                                      xformla = ~Population+CO2,
                                      data = did_master_cities_panel_no2_robustness
)
simple.attgt_retirement50k_uni_cov_robustness <- aggte(attgt_retirement50k_uni_cov_robustness, type = "simple")
dynamic.attgt_retirement50k_uni_cov_robustness <- aggte(attgt_retirement50k_uni_cov_robustness, type = "dynamic")

mod1_no2_retirement50k_var_robustness <- tidy(dynamic.attgt_retirement50k_var_robustness) %>% mutate(model = "50k from Retired Coal Unit (Mod1)")
mod1_no2_retirement50k_var_cov_robustness <- tidy(dynamic.attgt_retirement50k_var_cov_robustness) %>% mutate(model = "50k from Retired Coal Unit (Mod2)")
mod1_no2_retirement50k_uni_robustness <- tidy(dynamic.attgt_retirement50k_uni_robustness) %>% mutate(model = "50k from Retired Coal Unit (Mod3)")
mod1_no2_retirement50k_uni_cov_robustness <- tidy(dynamic.attgt_retirement50k_uni_cov_robustness) %>% mutate(model = "50k from Retired Coal Unit (Mod4)")

mod1_no2_retirement_retonly <- bind_rows(mod1_no2_retirement50k_var_robustness, mod1_no2_retirement50k_var_cov_robustness,
                                mod1_no2_retirement50k_uni_robustness, mod1_no2_retirement50k_uni_cov_robustness)

pdf(file = "si_figures/figSI23_no2_main_retiredonly.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot(mod1_no2_retirement_retonly, aes(x = as.factor(event.time), y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange(aes(color = model, group = model), position = position_dodge(width=0.3)) +
  xlab("Time to Coal Power Unit Retirement") +
  ylab("Coefficient") +
  ggtitle("Average Effect by Length of Exposure") +
  theme_classic() + scale_color_brewer(palette = "Set1")  + 
  geom_vline(xintercept = "0", linetype = "twodash", color = "grey50", size=0.5) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", size=0.3)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S12: Staggered differences-in-differences, event study (NO2 Rates) by capacity.                             #                                    # 
#--------------------------------------------------------------------------------------------------------------------#

attgt_retirement50k_high_var_robustness <- att_gt(yname = "Rates_NO2",
                                       tname = "year",
                                       idname = "ID",
                                       gname = "Retired.year50k_high",
                                       base_period = "varying",
                                       bstrap = T,
                                       clustervars = "ID",
                                       control_group = "nevertreated",
                                       #xformla = ~Population+CO2,
                                       data = did_master_cities_panel_no2_robustness
)
simple.attgt_retirement50k_high_var_robustness <- aggte(attgt_retirement50k_high_var_robustness, type = "simple")
dynamic.attgt_retirement50k_high_var_robustness <- aggte(attgt_retirement50k_high_var_robustness, type = "dynamic")
cohort.attgt_retirement50k_high_var_no2_rate_robustness <- aggte(attgt_retirement50k_high_var_robustness, type = "group")

attgt_retirement50k_low_var_robustness <- att_gt(yname = "Rates_NO2",
                                      tname = "year",
                                      idname = "ID",
                                      gname = "Retired.year50k_low",
                                      base_period = "varying",
                                      bstrap = T,
                                      clustervars = "ID",
                                      control_group = "nevertreated",
                                      #xformla = ~Population+CO2,
                                      data = did_master_cities_panel_no2_robustness
)
simple.attgt_retirement50k_low_var_robustness <- aggte(attgt_retirement50k_low_var_robustness, type = "simple")
dynamic.attgt_retirement50k_low_var_robustness <- aggte(attgt_retirement50k_low_var_robustness, type = "dynamic")
cohort.attgt_retirement50k_low_var_no2_rate_robustness <- aggte(attgt_retirement50k_low_var_robustness, type = "group")

mod1_no2_retirement50k_high_var_robustness <- tidy(dynamic.attgt_retirement50k_high_var_robustness) %>% mutate(model = "50k from Retired Coal Unit (High Capacity MW)")
mod1_no2_retirement50k_low_var_robustness <- tidy(dynamic.attgt_retirement50k_low_var_robustness) %>% mutate(model = "50k from Retired Coal Unit (Low Capacity MW)")

mod1_no2_retirement_capacity_robustness <- bind_rows(mod1_no2_retirement50k_high_var_robustness, mod1_no2_retirement50k_low_var_robustness)

pdf(file = "si_figures/figSI24_no2_capacity_retiredonly.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
ggplot(mod1_no2_retirement_capacity_robustness, aes(x = as.factor(event.time), y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange(aes(color = model, group = model), position = position_dodge(width=0.3)) +
  xlab("Time to Coal Power Unit Retirement") +
  ylab("Coefficient") +
  ggtitle("Average Effect by Length of Exposure") +
  theme_classic() + scale_color_brewer(palette = "Set1")  + 
  geom_vline(xintercept = "0", linetype = "twodash", color = "grey50", size=0.5) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", size=0.3)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S13: FEct for Mortality Rates from PM2.5.                                                                   # 
#--------------------------------------------------------------------------------------------------------------------#

out.fect1 <- fect(Rates_PM ~ treated_after_treat50k , data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1),
                  index = c("ID","year"), 
                  method = "mc", force = "two-way", se = TRUE, parallel = TRUE, nboots = 500)

pdf(file = "si_figures/figSI7_fect_pm_rates.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
plot(out.fect1, ylab = "Effect of Coal Retirement on PM2.5 Rates", 
     cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S14: FEct for Emissions of PM2.5.                                                                           # 
#--------------------------------------------------------------------------------------------------------------------#

out.fect3 <- fect(PM ~ treated_after_treat50k , data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1), 
                  index = c("ID","year"), 
                  method = "mc", force = "two-way", se = TRUE, parallel = TRUE, nboots = 500)

pdf(file = "si_figures/figSI9_fect_pm_emissions.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
plot(out.fect3, ylab = "Effect of Coal Retirement on PM2.5 Emissions", 
     cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S15: FEct for Mortality Rates from NO2.                                                                     # 
#--------------------------------------------------------------------------------------------------------------------#

out.fect2 <- fect(Rates_NO2 ~ treated_after_treat50k , data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1), 
                  index = c("ID","year"), 
                  method = "mc", force = "two-way", se = TRUE, parallel = TRUE, nboots = 500)

pdf(file = "si_figures/figSI8_fect_no2_rates.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
plot(out.fect2, ylab = "Effect of Coal Retirement on NO2 Rates", 
     cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S16: FEct for Emissions of NO2.                                                                             # 
#--------------------------------------------------------------------------------------------------------------------#

out.fect4 <- fect(NO2 ~ treated_after_treat50k , data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1), 
                  index = c("ID","year"), 
                  method = "mc", force = "two-way", se = TRUE, parallel = TRUE, nboots = 500)

pdf(file = "si_figures/figS10_fect_no2_emissions.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
plot(out.fect4, ylab = "Effect of Coal Retirement on NO2 Emissions", 
     cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S17: FEct for Mortality Rates from PM2.5.                                                                   # 
#--------------------------------------------------------------------------------------------------------------------#

out.fect1b <- fect(Rates_PM ~ treated_after_treat50k , data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1),
                   index = c("ID","year"), 
                   method = "fe", force = "two-way", se = TRUE, parallel = TRUE, nboots = 500)

pdf(file = "si_figures/figSI11_fect_pm_rates_fe.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
plot(out.fect1b, ylab = "Effect of Coal Retirement on PM2.5 Mortality Rates", 
     cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S18: FEct for Emissions of PM2.5.                                                                           # 
#--------------------------------------------------------------------------------------------------------------------#

out.fect3b <- fect(PM ~ treated_after_treat50k , data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1), 
                   index = c("ID","year"), 
                   method = "fe", force = "two-way", se = TRUE, parallel = TRUE, nboots = 500)

pdf(file = "si_figures/figSI13_fect_pm_emissions_fe.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
plot(out.fect3b, ylab = "Effect of Coal Retirement on PM2.5 Emissions", 
     cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S19: FEct for Mortality Rates from NO2                                                                      # 
#--------------------------------------------------------------------------------------------------------------------#

out.fect2b <- fect(Rates_NO2 ~ treated_after_treat50k , data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1), 
                   index = c("ID","year"), 
                   method = "fe", force = "two-way", se = TRUE, parallel = TRUE, nboots = 500)

pdf(file = "si_figures/figSI12_fect_no2_rates_fe.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
plot(out.fect2b, ylab = "Effect of Coal Retirement on NO2 Mortality Rates", 
     cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S20: FEct for Emissions of NO2.                                                                             # 
#--------------------------------------------------------------------------------------------------------------------#

out.fect4b <- fect(NO2 ~ treated_after_treat50k , data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1), 
                   index = c("ID","year"), 
                   method = "fe", force = "two-way", se = TRUE, parallel = TRUE, nboots = 500)

pdf(file = "si_figures/figSI14_fect_no2_emissions_fe.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
plot(out.fect4b, ylab = "Effect of Coal Retirement on NO2 Emissions", 
     cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8)
dev.off()

#--------------------------------------------------------------------------------------------------------------------#
# Figure S21: Sensitivity Analysis PM2.5 Mortality                                                                   # 
#--------------------------------------------------------------------------------------------------------------------#

master_cities_panel <- master_cities_panel %>%
  mutate(logco2 = log(CO2))

model1 <- lm(Rates_PM~treated_after_treat50k+logco2+as.factor(ID)+as.factor(year), data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))

sensitivity_pm_rates <- sensemakr(model = model1, 
                                  treatment = "treated_after_treat50k",
                                  benchmark_covariates = "logco2",
                                  kd = 1:4)
pdf(file = "si_figures/figSI4_sensitivity_pm_rates.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
plot(sensitivity_pm_rates)
dev.off()

ovb_minimal_reporting(sensitivity_pm_rates)

#--------------------------------------------------------------------------------------------------------------------#
# Figure S22: Sensitivity Analysis PM2.5 Emissions                                                                   # 
#--------------------------------------------------------------------------------------------------------------------#

model3 <- lm(PM~treated_after_treat50k+logco2+as.factor(ID)+as.factor(year), data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))

sensitivity_pm_emissions <- sensemakr(model = model3, 
                                      treatment = "treated_after_treat50k",
                                      benchmark_covariates = "logco2",
                                      kd = 1:4)

pdf(file = "si_figures/figSI5_sensitivity_pm_emissions.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
plot(sensitivity_pm_emissions)
dev.off()

ovb_minimal_reporting(sensitivity_pm_emissions)

#--------------------------------------------------------------------------------------------------------------------#
# Figure S23: Sensitivity Analysis NO2 Mortality                                                                     # 
#--------------------------------------------------------------------------------------------------------------------#

model2 <- lm(Rates_NO2~treated_after_treat50k+logco2+as.factor(ID)+as.factor(year), data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))

sensitivity_no2_rates <- sensemakr(model = model2, 
                                   treatment = "treated_after_treat50k",
                                   benchmark_covariates = "logco2",
                                   kd = 1:4)

pdf(file = "si_figures/figSI3_sensitivity_no2_rates.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
plot(sensitivity_no2_rates)
dev.off()

ovb_minimal_reporting(sensitivity_no2_rates)

#--------------------------------------------------------------------------------------------------------------------#
# Figure S24: Sensitivity Analysis NO2 Emissions                                                                     # 
#--------------------------------------------------------------------------------------------------------------------#

model4 <- lm(NO2~treated_after_treat50k+logco2+as.factor(ID)+as.factor(year), data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))

sensitivity_no2_emissions <- sensemakr(model = model4, 
                                       treatment = "treated_after_treat50k",
                                       benchmark_covariates = "logco2",
                                       kd = 1:4)

pdf(file = "si_figures/figSI6_sensitivity_no2_emissions.pdf",   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 10) # The height of the plot in inches
plot(sensitivity_no2_emissions)
dev.off()

ovb_minimal_reporting(sensitivity_no2_emissions)

###############################################
########## TABLES MAIN TEXT ANALYSIS ##########
###############################################

#All sample, country, city, year FE
mod1 <- feols(Rates_NO2~v2x_polyarchy+Population+log(CO2)|ID+country2+year, data = master_pollution_dem)
mod2 <- feols(Rates_PM~v2x_polyarchy+Population+log(CO2)|ID+country2+year, data = master_pollution_dem)

#Border version 1, country, city, year FE
mod3 <- feols(Rates_NO2~v2x_polyarchy+Population+log(CO2)|ID+country2+year, data = subset(master_pollution_dem, border_city_t1 == 1))
mod4 <- feols(Rates_PM~v2x_polyarchy+Population+log(CO2)|ID+country2+year, data = subset(master_pollution_dem, border_city_t1 == 1))

#Border version 2: country, city, year FE
mod5 <- feols(Rates_NO2~v2x_polyarchy+Population+log(CO2)|ID+country2+year, data = subset(master_pollution_dem, border_city_t2 == 1))
mod6 <- feols(Rates_PM~v2x_polyarchy+Population+log(CO2)|ID+country2+year, data = subset(master_pollution_dem, border_city_t2 == 1))

#Alternative version of democracy, border version 1
mod7 <- feols(Rates_NO2~political_regime+Population+log(CO2)|ID+country2+year, data = subset(master_pollution_dem, border_city_t1 == 1))
mod8 <- feols(Rates_PM~political_regime+Population+log(CO2)|ID+country2+year, data = subset(master_pollution_dem, border_city_t1 == 1))

#Alternative version of democracy, border version 1
mod9 <- feols(Rates_NO2~polity_iv_dem+Population+log(CO2)|ID+country2+year, data = subset(master_pollution_dem, border_city_t1 == 1))
mod10 <- feols(Rates_PM~polity_iv_dem+Population+log(CO2)|ID+country2+year, data = subset(master_pollution_dem, border_city_t1 == 1))

#Only city and year FE
mod11 <- feols(Rates_NO2~v2x_polyarchy+Population+log(CO2)|ID+year, data = subset(master_pollution_dem, border_city_t1 == 1))
mod12 <- feols(Rates_PM~v2x_polyarchy+Population+log(CO2)|ID+year, data = subset(master_pollution_dem, border_city_t1 == 1))

#--------------------------------------------------------------------------------------------------------------------#
# Table 1: Democracy and Air Pollution (Particulate Matter)                                                          # 
#--------------------------------------------------------------------------------------------------------------------#

#PM2.5
texreg(list(mod2, mod4, mod6, mod12, mod8, mod10))

#--------------------------------------------------------------------------------------------------------------------#
# Table 2: Democracy and Air Pollution (NO2)                                                                         # 
#--------------------------------------------------------------------------------------------------------------------#

#NO2
texreg(list(mod1, mod3, mod5, mod11, mod7, mod9))

#--------------------------------------------------------------------------------------------------------------------#
# Table 3: Climate Policy Commitment and Coal Phase Out                                                              # 
#--------------------------------------------------------------------------------------------------------------------#

model1a <- feols(log(tot_retired_capacity+1)~lag(tot_mitigation_policies)+NY.GNP.PCAP.CDK, data = subset(master_plants_panel, coal_country == 1))
model1b <- feols(log(tot_retired_plants+1)~lag(tot_mitigation_policies)+NY.GNP.PCAP.CDK, data = subset(master_plants_panel, coal_country == 1))
model1c <- feols(retirement_bin~lag(tot_mitigation_policies)+NY.GNP.PCAP.CDK, data = subset(master_plants_panel, coal_country == 1))

model1d <- feols(log(tot_retired_capacity+1)~lag(tot_mitigation_policies)+NY.GNP.PCAP.CDK|country2+year, data = subset(master_plants_panel, coal_country == 1))
model1e <- feols(log(tot_retired_plants+1)~lag(tot_mitigation_policies)+NY.GNP.PCAP.CDK|country2+year, data = subset(master_plants_panel, coal_country == 1))
model1f <- feols(retirement_bin~lag(tot_mitigation_policies)+NY.GNP.PCAP.CDK|country2+year, data = subset(master_plants_panel, coal_country == 1))

texreg(list(model1a, model1b, model1c, model1d, model1e, model1f))

##############################################
######## TABLES SUPPORTING INFORMATION #######
##############################################

master_cities_panel <- master_cities_panel %>%
  mutate(logco2 = log(CO2))

#--------------------------------------------------------------------------------------------------------------------#
# Table S2: Standard TWFE: Mortality Rates.                                                                          # 
#--------------------------------------------------------------------------------------------------------------------#

model1 <- feols(Rates_NO2~treated_after_treat50k+logco2|ID+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model2 <- feols(Rates_PM~treated_after_treat50k+logco2|ID+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))

texreg(list(model1, model2))

#--------------------------------------------------------------------------------------------------------------------#
# Table S3: Standard TWFE: Total Emissions.                                                                          # 
#--------------------------------------------------------------------------------------------------------------------#

model3 <- feols(NO2~treated_after_treat50k+logco2|ID+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model4 <- feols(PM~treated_after_treat50k+logco2|ID+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model5 <- feols(logco2~treated_after_treat50k|ID+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))

texreg(list(model3, model4, model5))

#--------------------------------------------------------------------------------------------------------------------#
# Table S4: Standard TWFE: Mortality Rates.                                                                          # 
#--------------------------------------------------------------------------------------------------------------------#

model6 <- feols(Rates_NO2~treated_after_treat50k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model7 <- feols(Rates_PM~treated_after_treat50k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))

texreg(list(model6, model7))

#--------------------------------------------------------------------------------------------------------------------#
# Table S5: Standard TWFE: Total Emissions                                                                           # 
#--------------------------------------------------------------------------------------------------------------------#

model8 <- feols(NO2~treated_after_treat50k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model9 <- feols(PM~treated_after_treat50k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model10 <- feols(logco2~treated_after_treat50k|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))

texreg(list(model8, model9, model10))

#--------------------------------------------------------------------------------------------------------------------#
# Table S6: Standard Two-Way Fixed Effects: Mortality Rates by Distance to the Closest Retired Unit                  #                                                      
#--------------------------------------------------------------------------------------------------------------------#

model11 <- feols(Rates_PM~treated_after_treat50k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model12 <- feols(Rates_PM~treated_after_treat100k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating100k == 1 | close_retired100k == 1))
model13 <- feols(Rates_PM~treated_after_treat150k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating150k == 1 | close_retired150k == 1))

texreg(list(model11, model12, model13))

#--------------------------------------------------------------------------------------------------------------------#
# Table S7: Standard Two-Way Fixed Effects: Mortality Rates by Distance to the Closest Retired Unit.                 #                                                       
#--------------------------------------------------------------------------------------------------------------------#

model14 <- feols(Rates_NO2~treated_after_treat50k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model15 <- feols(Rates_NO2~treated_after_treat100k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating100k == 1 | close_retired100k == 1))
model16 <- feols(Rates_NO2~treated_after_treat150k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating150k == 1 | close_retired150k == 1))

texreg(list(model14, model15, model16))

#--------------------------------------------------------------------------------------------------------------------#
# Table S8: Standard Two-Way Fixed Effects: Emissions by Distance to the Closest Retired Unit.                       #                                                  
#--------------------------------------------------------------------------------------------------------------------#

model17 <- feols(PM~treated_after_treat50k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model18 <- feols(PM~treated_after_treat100k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating100k == 1 | close_retired100k == 1))
model19 <- feols(PM~treated_after_treat150k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating150k == 1 | close_retired150k == 1))

texreg(list(model17, model18, model19))

#--------------------------------------------------------------------------------------------------------------------#
# Table S9: Standard Two-Way Fixed Effects: Emissions by Distance to the Closest Retired Unit.                       #                                               
#--------------------------------------------------------------------------------------------------------------------#

model20 <- feols(NO2~treated_after_treat50k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model21 <- feols(NO2~treated_after_treat100k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating100k == 1 | close_retired100k == 1))
model22 <- feols(NO2~treated_after_treat150k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating150k == 1 | close_retired150k == 1))

texreg(list(model20, model21, model22))

#--------------------------------------------------------------------------------------------------------------------#
# Table S10: Standard Two-Way Fixed Effects: Emissions by Distance to the Closest Retired Unit.                      #                                                     
#--------------------------------------------------------------------------------------------------------------------#

model23 <- feols(logco2~treated_after_treat50k|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model24 <- feols(logco2~treated_after_treat100k|ID+country2+year, data = subset(master_cities_panel, close_operating100k == 1 | close_retired100k == 1))
model25 <- feols(logco2~treated_after_treat150k|ID+country2+year, data = subset(master_cities_panel, close_operating150k == 1 | close_retired150k == 1))

texreg(list(model23, model24, model25))

#--------------------------------------------------------------------------------------------------------------------#
# Table S11: Standard Two-Way Fixed Effects: Mortality Rates by Capacity of the Closest Retired Coal Unit.           #                                                             
#--------------------------------------------------------------------------------------------------------------------#

model26 <- feols(Rates_PM~treated_after_treat50k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model27 <- feols(Rates_PM~treated_after_treat50k_high+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model28 <- feols(Rates_PM~treated_after_treat50k_low+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))

texreg(list(model26, model27, model28))

#--------------------------------------------------------------------------------------------------------------------#
# Table S12: Standard Two-Way Fixed Effects: Mortality Rates by Capacity of the Closest Retired Coal Unit.           #                                                                
#--------------------------------------------------------------------------------------------------------------------#

model29 <- feols(Rates_NO2~treated_after_treat50k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model30 <- feols(Rates_NO2~treated_after_treat50k_high+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model31 <- feols(Rates_NO2~treated_after_treat50k_low+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))

texreg(list(model29, model30, model31))

#--------------------------------------------------------------------------------------------------------------------#
# Table S13: Standard Two-Way Fixed Effects: Emissions by Capacity of the Closest Retired Coal Unit.                 #                                                        
#--------------------------------------------------------------------------------------------------------------------#

model32 <- feols(PM~treated_after_treat50k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model33 <- feols(PM~treated_after_treat50k_high+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model34 <- feols(PM~treated_after_treat50k_low+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))

texreg(list(model32, model33, model34))

#--------------------------------------------------------------------------------------------------------------------#
# Table S14: Standard Two-Way Fixed Effects: Emissions by Capacity of the Closest Retired Coal Unit.                 #                                                         
#--------------------------------------------------------------------------------------------------------------------#

model35 <- feols(NO2~treated_after_treat50k+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model36 <- feols(NO2~treated_after_treat50k_high+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model37 <- feols(NO2~treated_after_treat50k_low+logco2|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))

texreg(list(model35, model36, model37))

#--------------------------------------------------------------------------------------------------------------------#
# Table S15: Standard Two-Way Fixed Effects: Emissions by Capacity of the Closest Retired Coal Unit                  #                                                        
#--------------------------------------------------------------------------------------------------------------------#

model38 <- feols(logco2~treated_after_treat50k|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model39 <- feols(logco2~treated_after_treat50k_high|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))
model40 <- feols(logco2~treated_after_treat50k_low|ID+country2+year, data = subset(master_cities_panel, close_operating50k == 1 | close_retired50k == 1))

texreg(list(model38, model39, model40))

#--------------------------------------------------------------------------------------------------------------------#
# Table S16: Standard Two-Way Fixed Effects: Full Sample of All Cities.                                              #                            
#--------------------------------------------------------------------------------------------------------------------#

model41 <- feols(Rates_NO2~treated_after_treat50k+logco2|ID+country2+year, data = master_cities_panel)
model42 <- feols(Rates_PM~treated_after_treat50k+logco2|ID+country2+year, data = master_cities_panel)
model43 <- feols(PM~treated_after_treat50k+logco2|ID+country2+year, data = master_cities_panel)
model44 <- feols(NO2~treated_after_treat50k+logco2|ID+country2+year, data = master_cities_panel)

texreg(list(model41, model42, model43, model44))

#--------------------------------------------------------------------------------------------------------------------#
# Table S17: Standard Two-Way Fixed Effects: Full Sample of All Cities.                                              #                          
#--------------------------------------------------------------------------------------------------------------------#

model45 <- feols(Rates_NO2~treated_after_treat50k_high+logco2|ID+country2+year, data = master_cities_panel)
model46 <- feols(Rates_PM~treated_after_treat50k_high+logco2|ID+country2+year, data = master_cities_panel)
model47 <- feols(PM~treated_after_treat50k_high+logco2|ID+country2+year, data = master_cities_panel)
model48 <- feols(NO2~treated_after_treat50k_high+logco2|ID+country2+year, data = master_cities_panel)

texreg(list(model45, model46, model47, model48))
